Assignment7 Weather
Yunhong Feng
2024-02-28
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.5.0 ✔ purrr 0.3.4
## ✔ tibble 3.2.1 ✔ dplyr 1.1.4
## ✔ tidyr 1.1.4 ✔ stringr 1.4.0
## ✔ readr 2.1.1 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(dplyr)
library(ggplot2)
library(viridis)## Loading required package: viridisLitelibrary(readr)
library(sf)## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(maps)##
## Attaching package: 'maps'
##
## The following object is masked from 'package:viridis':
##
## unemp
##
## The following object is masked from 'package:purrr':
##
## maplibrary(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.library(leaflet)
library(RColorBrewer)
library(leaflet.extras) data<-read.csv("storms.csv")Task 1a) State Choropleth Maps
# Get total damage
data <- data %>%
mutate(TOTAL_DAMAGE = DAMAGE_PROPERTY_USD + DAMAGE_CROPS_USD)
# Standardize the format of STATE_FIPS for both data
damage_data_state <- data %>%
group_by(STATE_FIPS) %>%
summarise(TOTAL_DAMAGE = sum(TOTAL_DAMAGE, na.rm = TRUE))
damage_data_state <- damage_data_state %>% mutate(STATE_FIPS = sprintf("%02d", as.numeric(STATE_FIPS)))
# Get state level data and join with the weather data
# Since the number of TOTAL_DAMAGE is too large, I log it to make the scale more readable
states <- states(cb = TRUE, year = 2020, class = "sf")##
|
| | 0%
|
|= | 1%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|====== | 9%
|
|======= | 10%
|
|======== | 11%
|
|========= | 12%
|
|=========== | 16%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 22%
|
|================= | 24%
|
|================= | 25%
|
|==================== | 29%
|
|========================= | 35%
|
|========================== | 37%
|
|============================ | 40%
|
|================================ | 45%
|
|================================== | 48%
|
|==================================== | 51%
|
|===================================== | 52%
|
|====================================== | 55%
|
|========================================= | 58%
|
|========================================== | 60%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 66%
|
|================================================== | 71%
|
|==================================================== | 74%
|
|===================================================== | 76%
|
|====================================================== | 78%
|
|======================================================= | 79%
|
|========================================================= | 81%
|
|========================================================== | 82%
|
|============================================================= | 87%
|
|============================================================== | 88%
|
|=============================================================== | 90%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|==================================================================== | 97%
|
|======================================================================| 99%
|
|======================================================================| 100%map_data1a <- left_join(states, damage_data_state, by = c("GEOID" = "STATE_FIPS")) %>%
filter(!is.na(TOTAL_DAMAGE))
map_data1a$LOG_TOTAL_DAMAGE <- log(map_data1a$TOTAL_DAMAGE)
# Plot the map
ggplot(map_data1a) +
geom_sf(aes(fill = LOG_TOTAL_DAMAGE), color = "white") +
scale_fill_viridis_c(
name = "Total Damage (log10)") +
labs(title = "Total Damage from Storms by State",
subtitle = "United States",
caption = "Source: Storm Data")+
geom_sf_text(aes(label = STUSPS), color = "black", size = 3, check_overlap = TRUE) + # State code
theme_minimal() +
theme(legend.position = "bottom") +
coord_sf(xlim = c(-125, -66), ylim = c(24, 49), expand = FALSE) # Center and Zoom in US## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
Task 1b) County Choropleth Maps
# Repeat the process to create county level map
damage_data_county <- data %>%
group_by(CZ_FIPS) %>%
summarise(TOTAL_DAMAGE = sum(TOTAL_DAMAGE, na.rm = TRUE))
damage_data_county <- damage_data_county %>% mutate(CZ_FIPS = sprintf("%03d", as.numeric(CZ_FIPS)))
counties <- counties(cb = TRUE, year = 2020, class = "sf")##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 3%
|
|=== | 4%
|
|============ | 17%
|
|==================== | 29%
|
|===================== | 30%
|
|============================== | 43%
|
|================================== | 49%
|
|=========================================== | 62%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|======================================================================| 99%
|
|======================================================================| 100%map_data1b <- left_join(counties, damage_data_county, by = c("COUNTYFP" = "CZ_FIPS")) %>%
filter(!is.na(TOTAL_DAMAGE))
map_data1b$LOG_TOTAL_DAMAGE <- log(map_data1b$TOTAL_DAMAGE)
ggplot(map_data1b) +
geom_sf(aes(fill = LOG_TOTAL_DAMAGE), color = "white") +
scale_fill_viridis_c(name = "Total Damage (log10)") +
labs(title = "Total Damage from Storms by County",
subtitle = "United States",
caption = "Source: Storm Data") +
theme_minimal() +
theme(legend.position = "bottom") +
coord_sf(xlim = c(-125, -66), ylim = c(24, 49), expand = FALSE)
Task 1c) Density Map
# Data cleaning
Injuries_data <- data %>%
filter(!is.na(BEGIN_LAT) & !is.na(BEGIN_LON))%>%
group_by(BEGIN_LAT, BEGIN_LON) %>%
summarise(Total_Injuries = sum(INJURIES_DIRECT + INJURIES_INDIRECT, na.rm = TRUE)) %>%
ungroup() %>% filter(Total_Injuries > 0) # Only focus on the events that causes injury, make the graph more useful## `summarise()` has grouped output by 'BEGIN_LAT'. You can override using the
## `.groups` argument.# Get US map
storm_points <- st_as_sf(Injuries_data, coords = c("BEGIN_LON", "BEGIN_LAT"), crs = 4326)
us_map <- st_as_sf(maps::map("usa", plot = FALSE, fill = TRUE), crs = 4326)
# Plot the dots that shows the density on the US map
ggplot()+
geom_sf(data = states, fill = "lightgrey", color = "black", size = 0.25) +
geom_sf(data = storm_points, aes(size = Total_Injuries, color = Total_Injuries), alpha = 0.5) +
scale_color_viridis(option = "C", direction = -1) + # Use viridis for a nice color scale
theme_minimal() +
labs(title = "Density of Severe Weather Events by Injuries",
subtitle = "United States",
color = "Injuries",
size = "Injuries",
caption = "Source: Storm Data") +
theme(legend.position = "bottom") +
coord_sf(xlim = c(-130, -65), ylim = c(25, 50)) # Center and Zoom in US
Both the Choropleth Maps and the Density Map have pros and cons in their
visualization. The County Choropleth Map excellently visualizes how the
economic impacts of storms distribute across larger geographic areas,
highlighting regions that are more financially affected. Whereas the
density map can more precisely indicate the exact locations of severe
weather events, providing a clear picture of where the deadliest storms
occur. From my perspective the point-based map may be more effective due
to its direct representation of fatalities, offering a stark
visualization of the human cost of storms.
Task 2a) Interactive Map of Severe Weather Events
# Select useful data info
leaflet_data<- data %>%
mutate(Total_Deaths = DEATHS_DIRECT + DEATHS_INDIRECT,
Total_Injuries = INJURIES_DIRECT + INJURIES_INDIRECT) %>%
filter(Total_Deaths > 0) %>%
select(BEGIN_LAT, BEGIN_LON, BEGIN_DATE_TIME, END_DATE_TIME, EPISODE_ID, EVENT_ID, Total_Deaths, Total_Injuries, EVENT_TYPE)
# Plot the graph with popup
leaflet(data = leaflet_data) %>%
addMarkers(lng = ~BEGIN_LON, lat = ~BEGIN_LAT, popup = ~paste(
"Begin Date: ", BEGIN_DATE_TIME,
"<br>End Date: ", END_DATE_TIME,
"<br>Episode ID: ", EPISODE_ID,
"<br>Event ID: ", EVENT_ID,
"<br>Event Type: ", EVENT_TYPE,
"<br>Total Deaths: ", Total_Deaths,
"<br>Total Injuries: ", Total_Injuries)) %>%
addTiles(urlTemplate = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png") %>%
setView(lng = -98.5795, lat = 39.8283, zoom = 4) # Center the map over the US## Warning in validateCoords(lng, lat, funcName): Data contains 2335 rows with
## either missing or invalid lat/lon values and will be ignoredTask 2b) Color by Type of Weather Event
# Determine what are the types of event in the data
table(leaflet_data$EVENT_TYPE)##
## Astronomical Low Tide Avalanche Blizzard
## 1 71 27
## Coastal Flood Cold/Wind Chill Debris Flow
## 1 125 10
## Dense Fog Dust Storm Excessive Heat
## 50 11 265
## Extreme Cold/Wind Chill Flash Flood Flood
## 76 249 129
## Freezing Fog Frost/Freeze Hail
## 3 2 7
## Heat Heavy Rain Heavy Snow
## 593 29 35
## High Surf High Wind Hurricane
## 59 45 59
## Ice Storm Lake-Effect Snow Lakeshore Flood
## 11 6 1
## Lightning Marine Dense Fog Marine Strong Wind
## 106 4 9
## Marine Thunderstorm Wind Marine Tropical Storm Rip Current
## 8 3 380
## Sleet Sneakerwave Storm Surge/Tide
## 2 9 11
## Strong Wind Thunderstorm Wind Tornado
## 75 189 109
## Tropical Depression Tropical Storm Waterspout
## 3 50 1
## Wildfire Winter Storm Winter Weather
## 76 70 211# Re - categorize the data into simpler form
leaflet_data <- leaflet_data %>%
mutate(EVENT_TYPE2 = case_when(
EVENT_TYPE %in% c("Heat", "Cold/Wind Chill", "Excessive Heat", "Extreme Cold/Wind Chill", "Frost/Freeze") ~ "Extreme Temperature",
EVENT_TYPE %in% c("Flash Flood", "Flood", "Coastal Flood", "Lakeshore Flood", "Storm Surge/Tide") ~ "Flood",
EVENT_TYPE %in% c("High Wind", "Hurricane", "Tropical Storm", "Marine Strong Wind", "Strong Wind", "Thunderstorm Wind", "Tornado") ~ "Wind",
EVENT_TYPE %in% c("Blizzard", "Heavy Snow", "Ice Storm", "Lake-Effect Snow", "Winter Storm", "Winter Weather", "Sleet", "Freezing Fog") ~ "Winter Weather",
TRUE ~ "Other")) # Other unclassified events
# Assign each category a color
colors <- c("Extreme Temperature" = "red",
"Flood" = "blue",
"Wind" = "green",
"Winter Weather" = "cyan",
"Other" = "grey")
leaflet_data$Color <- sapply(leaflet_data$EVENT_TYPE2, function(type) colors[type])
# Plot the map with the color fill
leaflet(leaflet_data) %>%
addTiles() %>%
addCircleMarkers(lng = ~BEGIN_LON, lat = ~BEGIN_LAT,
color = ~Color, # Use the color column
popup = ~paste(
"Begin Date: ", BEGIN_DATE_TIME,
"<br>End Date: ", END_DATE_TIME,
"<br>Episode ID: ", EPISODE_ID,
"<br>Event ID: ", EVENT_ID,
"<br>Event Type: ", EVENT_TYPE,
"<br>Total Deaths: ", Total_Deaths,
"<br>Total Injuries: ", Total_Injuries),
radius = 6, fillOpacity = 0.8) %>%
setView(lng = -98.5795, lat = 39.8283, zoom = 4) %>%
addLegend(position = "bottomright", colors = colors, labels = names(colors), title = "Event Type")## Warning in validateCoords(lng, lat, funcName): Data contains 2335 rows with
## either missing or invalid lat/lon values and will be ignoredTask 2c) Cluster
# Add cluster code to the map
leaflet(data = leaflet_data) %>%
addTiles() %>%
addCircleMarkers(
lng = ~BEGIN_LON, lat = ~BEGIN_LAT,
popup = ~paste(
"Begin Date: ", BEGIN_DATE_TIME,
"<br>End Date: ", END_DATE_TIME,
"<br>Episode ID: ", EPISODE_ID,
"<br>Event ID: ", EVENT_ID,
"<br>Event Type: ", EVENT_TYPE,
"<br>Total Deaths: ", Total_Deaths,
"<br>Total Injuries: ", Total_Injuries),
clusterOptions = markerClusterOptions(),
group = 'markers') %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
setView(lng = -98.5795, lat = 39.8283, zoom = 4) ## Warning in validateCoords(lng, lat, funcName): Data contains 2335 rows with
## either missing or invalid lat/lon values and will be ignored